On two continuous models for the dynamics of sandpile surfaces 
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We consider a modified BCRE model for pile surface dy- 
namics and show that in the long-scale limit this model con- 
verges to a quasistationary model of pile growth in the form 
of an evolutionary variational inequality. 
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I. INTRODUCTION 

Recent much interest to the physics of granular me- 
dia was, in particular, stimulated by two salient fea- 
tures of the granular state: multiplicity of metastable 
pile shapes and occurrence of avalanches upon pile sur- 
faces. It has been realized that, to account for metasta- 
bility, the model of pile surface dynamics should not be 
written as an evolutionary equation for the pile surface 
alone. An additional unknown characterizing the flow of 
grains down the pile surface is useful because such flows 
are not uniquely determined by the external source and 
local free surface topography. 

A large spatio-temporal scale pile growth model involv- 
ing two coupled dependent variables and able to account 
for metastability has been proposed in 0J^]. This model 
neglects avalanches as small fluctuations of the pile sur- 
face and describes the evolving mean surface of a pile 
that grows on an arbitrary support under a given dis- 
tributed source of bulk material. The model permits an 
equivalent formulation as an evolutionary variational or 
quasivariational inequality; such a formulation simplifies 
significantly both the mathematical study of the prob- 
lem H and its numerical solution |l],f|,^) ■ As it has been 
shown in the shapes of real piles on flat open plat- 
forms [|| are described by the analytical solutions of this 
inequality. A modification of the model, able to account 
for avalanches as almost instantaneous slides, is also dis- 
cussed in Q ; according to observations made in the same 
work, such a slide may, indeed, be a possible avalanche 
scenario (see also 0). 

Independently and using different arguments, the same 
pile growth model in the form of a variational inequality 
has been derived by Aronsson, Evans, and Wu S. In 
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||, Evans et al. studied its discontinuous solutions cor- 
responding to avalanches; in [fl0[ Evans and Rezakhanlou 
showed that the cellular automata models of sandpiles, 
presented as intuitively attractive examples in almost all 
works on self-organized criticality flTl[ , converge in a con- 
tinuous limit to a similar variational inequality with an 
anisotropy inherited from the cellular structure of these 
crude models. 

A different continuous model, also involving two cou- 
pled dependent variables and describing the granular sur- 
face flow and pile surface dynamics, has been proposed 
by Bouchaud, Cates, Ravi Prakash, and Edwards (the 
BCRE model) fljjlj. Although the choice of the ba- 
sic variables in this model is equivalent to that in 0,0, 
the model is written for the free surfaces only slightly 
deviating from the critical slope and employs different 
phenomenological constitutive relations. The emphasis 
is put onto the simulation of fast processes, like ampli- 
fication and distinction of rolling grains population dur- 
ing an avalanche. The BCRE model has been simplified 
by de Gennes applied to various one-dimensional 

surface flow problems (see, e.g. , fli^ , |l6|| ), and modified 
for thick surface granular flows Jl7|. Further exact solu- 
tions to simplified BCRE equations can be constructed 
by the methods proposed in | 18 | . Using the BCRE 
model, Bouchaud and Cates |19| explained another type 
of avalanches (in a thin granular layer on an inclined 
plane, see pQ]). 

Our aim here is to investigate a relation between the 
two models mentioned above. After reminding briefly 
the variational and BCRE models, we propose a full- 
dimensional generalization of the latter, originally for- 
mulated by BCRE in the one-dimensional case. To do 
this, we modify and extend the constitutive relations de- 
termining the surface flow velocity and the rolling-to- 
immobilized-state transition rate: BCRE's assumption 
that the slope is everywhere almost critical is too restric- 
tive for our purpose. Rescaling the variables, we show 
that the modified BCRE model contains a small param- 
eter, the ratio of a characteristic rolling grains layer thick- 
ness to the pile size, and hence may often be simplified 
by employing a quasistationary equation for the rolling 
grains layer. The issue of scaling turns out to be very 
important in description of pile growth: another dimen- 
sionless parameter in the model thus obtained is the ratio 
of a typical rolling grain path length to the pile size. For 
large piles, this coefficient is also small and we show that 
in the long-scale limit the modified BCRE model tends 
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to the variational model |],0]. For small piles, the cor- 
responding term can be significant. These results make 
clear why differ the shapes of small and large piles and, 
correspondingly, why different models should be used to 
simulate, say, formation of large sand dunes and small 
Aeolian ripples. 



II. VARIATIONAL MODEL OF PILE GROWTH 

Let a cohesionless granular material having an angle of 
repose a r be tipped out onto a given rough rigid surface 
y = h (x), where x — (xi, x 2 ) G M 2 . We want to find the 
shape of a pile thus generated. 

The real process of pile growth is often intermittent: 
discharged granular material not only flows continuously 
over the pile slopes but is also able to build up and then to 
pour suddenly down the slope in an avalanche. However, 
the avalanches usually involve only a small amount of 
particles in a pile and cause small fluctuations of the pile 
free surface. The model 0,0] neglects these fluctuations 
and is a model for the mean surface evolution. Whether 
the pile evolution is governed by a continuous surface 
flow or results from many small avalanches, the surface 
flow is typically confined to a thin boundary layer which 
is distinctly separated from the motionless bulk . 

Let us assume for simplicity that the support surface 
has no steep slopes, i.e., 

\Vho\ < k, 

where k = tana r (see [|l]-[| for the general case). Assum- 
ing the bulk density of material in a pile is constant we 
can write the conservation law as 

d t h + V • q = w, 

where h(x,t) is the free surface, q(x,t) is the horizontal 
projection of the flux of rolling particles, and w(x,t) - 
the source intensity. We neglect the inertia and suppose 
that surface flow is directed towards the steepest descent, 



where 



-rriVh, 



m(x, t) > 



(1) 



is an unknown scalar function. The conservation law 
takes now the form 



d t h - V • (mVh) = w. 



(2) 



It is assumed in this model that the surface slope angle 
cannot exceed the angle of repose, 



|V7i| < k, 



(3) 



and that no pouring occurs over the parts of the pile 
surface which arc inclined less: 



\Vh(x,t)\ < k =^> m(x,t) = 0. (4) 
To complete the model we have to specify the initial, 



h\t=o — ho, 



(5) 



and a boundary condition. Let the granular material be 
allowed to leave the system freely through part Ti of 
the boundary of domain C ffi 2 , and the other part of 
the boundary, T2, presents an impermeable wall. The 
boundary conditions are then, respectively, 



h\rt = h \ Tl , md n h\i 



0. 



(G) 



The model (0)-(|) contains two coupled unknowns, the 
free surface h and an auxiliary function m determining 
the rolling grains flux magnitude. Conditions ([!]), (||), 
and Jj4j) define to as a multivalued function of \Vh\, see 
Fig. 



k 

FIG. 1. Multivalued constitutive relation m = m(\Vh\). 

The problem ([!])- @ may be considered an anoma- 
lous diffusion problem and solved by approximating this 
highly nonlinear multivalued relation. However, a bet- 
ter way to solve this problem is based on its following 
reformulation in the form of an evolutionary variational 
inequality (see |2^] and for variational inequalities 
in mechanics and physics and their numerical solution, 
respectively). 

Let us define the set K of possible surfaces as 

K = { v {x) I \V V \<k, <p\ Tl =h \ ri } 

and the scalar product of two functions as (</), ip) = 
Jq (j)ip dx. We can now consider the following problem 
(variational inequality): 

Find h(x,t) such that h G K for all t > 0, 

(d t h -w,<p-h)>0 for all ip G K, (7) 
and h\ t= o = hy. 

Theorem. Function h(x, t) is a solution of the vari- 
ational inequality (^) if and only if there exists m(x,t) 
such that the pair {h,m} is a solution to (Q)-(|). 

The outline of the proof is given in |2j (see for math- 
ematical details and a proof of existence of a unique 
solution to the variational inequality (^)). It has been 
also shown that the surface flux magnitude m(x, t) is, in 
this model, a Lagrange multiplier related to the point- 
wise constraint (|3|). The values of such multipliers are 
not uniquely determined by the local conditions, which 
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is the " mathematical explanation" of long-range interac- 
tions typical of extended dissipative systems in a critical 
(marginally stable) state, see ]24| . 

The model (0)-ra) or, equivalently, (0), has simple ana- 
lytical solutions HjV such as the conical pile growing under 
a point-like source or the piles on flat open platforms de- 
scribed in H . Numerical solutions (lJQ| also demonstrate 
simple geometrical structures that agree with one's sand- 
box memories. Although this model is much simplified 
in many respects, it allows for the multiplicity of possi- 
ble pile shapes. The avalanches may be introduced into 
the model as solution discontinuities (in time) triggered 
by sudden changes of the admissible set K and are 
instantaneous events. On the time scale of a slow pile 
growth the life of an avalanche is, indeed, very short. 



III. MODIFIED BCRE EQUATIONS 



The BCRE equations [|12|,[13| involve two coupled vari- 
ables: the pile height, y = h(x,t), and the effective 
thickness (density) of the rolling grains layer, R(x, t) 
(R(x, t)dfl is the volume that the material, currently 
rolling above the area dfi, would occupy in the pile). 
The model has been formulated for a two-dimensional 
pile (let 1 ); free surface slope deviations from the crit- 
ical angle were assumed small. Original BCRE equa- 
tions included diffusion terms to account for a non- 
locality of grains dislodgement and for fluctuations of 
rolling grains velocity. Although diffusion plays a crucial 
role in BCRE's scenario of avalanches |l^,[l^,|l{| , these 
terms were regularly omitted by other researchers who ei- 
ther assumed that in their problems diffusion is insignif- 
icant and simplified the model, or proposed a different 
avalanche scenario (see, e.g., jl4],[l6 18 



Below, we 

also omit the diffusion terms at first but introduce small 
diffusion at a later stage as a means for model regular- 
ization in transition to a large-scale limit. 

Simplified BCRE equations may be written as follows: 



d t h = T[h,R], d t R + d x {vR) 



T[h,R]. 



Here the term T[h, R] accounts for the conversion of 
rolling grains into immobilized grains and vice versa, v 
is the horizontal projection of rolling grains velocity, and 
w(x,t) is the source intensity (we assume that the tipped 
grains do not stick to the pile surface but join the rolling 
grains first). 

Limiting their consideration to the slopes that are close 
to critical, BCRE assumed constant downslope drift ve- 
locity v. The surface flux magnitude, q — vR, is thus 
determined solely by the rolling layer thickness R. Since 
in the previous model q = mk for the critical slopes, m 
and R play similar roles and the two choices of basic vari- 
ables, {h,m} and {h, R}, are essentially equivalent. The 
exchange term T in BCRE model is linearized in a vicin- 
ity of the critical angle a r and is proportional (for thin 



surface flows) to R: T[h, R] = r yR(a r — 9), where 9(x,t) 
is the surface slope angle and 7 is a coefficient. 

For a three-dimensional pile (1 e I 2 ), the model equa- 
tions are similar, 



d t h = T[h,R}, d t R + V • (vR) = w -T[h,R] 



(8) 



but the constitutive relations determining v and, prob- 
ably, T[h, R] should be modified; here we will follow 
p6[ (see also p7[|). We assume that the rolling parti- 
cles drift towards the steepest descent of the free surface 
with a mean velocity v depending on the slope angle (the 
steeper the slope, the higher is the velocity). On their 
way downslope, these particles may be trapped and ab- 
sorbed into the motionless bulk (the steeper the slope, 
the lower is the trapping rate T). If the surface is hor- 
izontal, the mean flow velocity is zero and the trapping 
rate is maximal, for 9 — a r the rolling particles follow 
without trapping. Below, we will not consider the over- 
critical slopes and assume also that the trapping rate is 
proportional to the amount of rolling grains R. 

At least partially, this simplified picture can be justi- 
fied by recent experimental, theoretical, and numerical 
studies on the motion of a spherical particle on a rough 
inclined plane [^8|,^9| . For the relevant region of slope an- 
gles 9, the energy dissipation due to the multiple shocks 
experienced by a moving particle is equivalent to the ac- 
tion of a viscous friction force pj| . Because of that such 
particles reach a constant mean velocity proportional to 
sin 9. Sometimes, however, the particles are suddenly 
trapped in a well and completely lose their momentum 
in the direction of motion (^9|. Of course, conditions 
in the collective flow of grains over the pile surface are 
somewhat different. In particular, the flow velocity may 
depend on the thickness of rolling grains layer |30f| and 
the exchange rate is not exactly proportional to R [ p"7| . 
Various improved dependencies can be incorporated into 
the model. The limiting behavior of the modified BCRE 
model is, however, robust and does not depend on de- 
tails. For clarity of presentation we will consider the 
long-scale limit of a thin-flow model with the simplest 
phenomenological relations determining the flow velocity 
and rolling-to-immobilized-state transition. 

Since the mean velocity of surface flow is proportional 
to sin 9 psj ], its horizontal projection v is proportional 
to sin 9 cos 9 = tan 9/(1 +tan 2 #). Postulating that the 
flow is in the steepest descent direction, we obtain v = 
— /xV/i/(l + |V/i| 2 ), where \i is a coefficient. Simplifying 
this relation we assume 



-fiVh. 



(9) 



The exchange rate T should not depend on the slope ori- 
entation and we assume it to be a smooth decreasing 
function of \Vh\ 2 that becomes zero for critical slopes. 
Assuming T is proportional to R (thin flows) we arrive 
at 



T[h, R] = 7.R 1 



|V/i|< 



(10) 
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as the simplest constitutive relation p6| . We will now de- 
rive a dimensionless formulation for the modified BCRE 
model (§)-©• 

The parameters in this model have the following di- 
mensions: [7] = T _1 and = CI . Let us denote 
by W the characteristic intensity of the external source; 
[w] = £T _1 . The three length scales characterizing the 
pile surface dynamics and surface granular flow may be 
defined as follows: 

• typical thickness of the rolling grains layer, Lp = 
W/r, 

• mean path of a rolling particle before it is trapped 
strongly depends on the slope steepness but, for a 
fixed subcritical slope, is proportional to the ratio 
Lp = fi/j characterizing the competition between 
rolling and trapping; 

• the pile size L. 

The time T — L/w needed for a source with given 
intensity w to produce a pile of size L may be used as a 
long time scale. Rescaling the variables, 

x = —x, h! = —h, R,' = — — R, w' = —w, t' = —t, 
L L Lp w 1 

we arrive at the following dimensionless formulation: 

9th = T[h,R], (11) 



L 



-d t R 



L, 



V • (RVh) 



T[h,R], 



T[h, R] = R[l 



Nh\ [ 



(12) 



(13) 



Typically, Lp <C Lp < L. The first coefficient in jl2] ) 
is very small, so it may often be possible to omit the 
corresponding term and use a quasistationary equation 
for the rolling layer. Such an approach has already been 
employed in simulation of the dynamics of sand ripples, 
see p6fl . The second coefficient, Lp/L, may be significant 
for small piles, like sand ripples, but becomes small too 
for large piles. Further simplification of the model is then 
appropriate. 



IV. THE LONG-SCALE LIMIT OF BCRE MODEL 

Let us denote v = Lp/L and study the v — > behavior 
of the model (|ll|)-([l^). This limit corresponds to the 
case of large piles (L ^> Lp). We want to show that 
in this limit the pile shape evolution is described by the 
variational inequality ([?]) which remains invariant under 
the rescaling employed. 

Physically, the situation is clear: although the model 
(|ll|)-(0) permits grains to roll down upon any inclined 



slope, the rolling particles are quickly stopped and their 
paths are short comparing to the pile size for all except 
the almost critical slopes. This is essentially what is 
assumed in the model which permits rolling upon 
the critical slopes only. Mathematically, the situation is 
somewhat more complicated. 

Since Lp <C Lp, we assume Lp/L is o(y) and set 
Lp/L = v\{v), where A tends to zero as v — > 0. Let us in- 
troduce a new variable, m = vR, define ip(u) = l — u 2 /k 2 , 
and rewrite the model ([llf)-([i3D as 



d t h=—^ ^, A9 t m-V • (mVTi) =w — -. 

v v 

For any v > this system consists of two coupled hy- 
perbolic equations. The second equation, which can be 
regarded as an equation for m, contains in its main part 
the coefficient V/i which may be discontinuous. The the- 
ory for such equations is complicated and not well devel- 
oped. To circumvent the difficulty, we add small diffusion 
to both equations and consider the regularized model 



dth 



rmp(\Vh\) 



e h Ah, 



(14) 



Tnib { V/z- ) 

\d t m - V • (mVh) = w — - + e m Am, (15) 

where the positive coefficients £h{v) and e m (v) vanish as 
v tends to zero. It should be noted that, although small 
diffusion may be physically meaningful and has been in- 
cluded into the original BCRE formulation |l2|,[l3|, here 
we introduce it merely as a parabolic regularization of hy- 
perbolic equations convenient for analyzing the model's 
behavior at v — > Q. 

We assume the same initial and boundary conditions, 
(||) and (|^) correspondingly, for the function h. The non- 
negative values of m both in Q at t = and on the bound- 
ary of this domain for t > may be chosen arbitrary: 
these initial and boundary conditions result only in the 
appearance of boundary layers in the solution for any fi- 
nite v > and are lost in the v — > limit. Rigorously, 
convergence of the problem (jl^)-(K) to the variational 
inequality (Q) is proved elsewhere Here we present 
a simplified scheme of the proof and avoid technicalities. 

The main step is, as usual, to obtain uniform in v > 
a priori estimates on the solutions of the equations 
(III)" (HI)- First, taking the gradient and multiplying by 
V/i, we derive from (|l4j) a parabolic partial differential 
equation for |V/i| 2 . Since ip(k) — and |V7io| < k, we 
are able, using the maximum principle for this equation, 
to show that for v > 

|V7i| < k for all (x, t); h is uniformly bounded. (16) 

Second, using the non-negativeness of the source func- 
tion w(x, t) and applying the maximum principle to the 
equation (|l5|), we deduce that for each v > 
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m > for all (x, t). 



(17) 



Applying the estimates (|l6|), ( |l7| ) to the equation (14) 
we obtain 



mil>(\Vh\) = 0(f) 



(18) 



Sending v to zero in (|l6|), (|l7|), ( |T8| ) we establish the 
fulfillment in this limit of the conditions rtl]) , (||), and (|J). 
Finally, adding the equations (jl4|) and ( Jl5| ) we obtain 

Xdtm + dth — V • (mV/i) = to + ShAh + e m Am. 

Since A, £/,, and e rn vanish as v — > 0, we can show that the 
corresponding limits of h and m satisfy also the balance 
equation (jsj) in some weak (integral) sense. This com- 
pletes the proof, because the model ((l])-(||) is equivalent 
to the variational inequality (0). 

To illustrate this result we will now compare solutions 
of the BCRE-type model (|l|)-((l|), solutions of the vari- 
ational inequality ([?]), and real shapes of small and large 
piles. Let us consider the simplest situation: a pile grow- 
ing under a point source on an infinite horizontal support 
Iiq = 0. Although in this case the piles are known to 
be almost perfect cones, sometimes one can notice J32] 
curved tails near the bottom of a small pile (Fig. |2pJ. 
As the pile becomes larger, the tail remains of only, say, 
tens of grain diameters long, so the tail of a large pile is 
difficult to see (Fig. ||b). 



(a) 



FIG. 2. Piles growing below a point-like source, (a) Dig- 
itized image of a small polenta pile, Alonso and Herrmann 
[32]. Different gray levels show the pile at different stages of 
growth, the curved tails near the pile bottom remain short; 
(b) Large conical pile, photo by Slater [33]. 



The modified BCRE model (|lJ)-lL5j) describes this sit- 
uation quite satisfactory, see Fig. |3J Although the tails 
of small piles [y = 0.2) are clearly seen, tails of the larger 
piles (v — 0.01) is difficult to detect. We see also that 
piles, predicted by the BCRE model with small v and A, 
are very close to the growing cone, exact analytical so- 
lution of the variational inequality (^). It may be noted 
that for small values of v and A the model equations are 
stiff and their numerical solution becomes difficult. Thus, 
even using an implicit finite-difference approximation of 
(|I^)-(|l5|), we needed 10 5 time steps in the latter example. 




(b) 





FIG. 3. Piles growing below a point-like source, numerical 
solution of (Q-(|l|) with e h = e m = 0, mj i=0 = 0, A = O.lv, 
and k = 1. 

(a) Small piles, v = 0.2; (b) Solid line - large piles, v = 0.01; 
dashed line - analytical solution of the variational inequality 

(I). 




V. CONCLUSION 

We considered two different continuous models for the 
pile surface dynamics: the BCRE model Jl2],[l3| and the 
variational model j[],U ■ Both models are written for two 
coupled dependent variables and are able, in principle, 
to account for multiplicity of metastable pile shapes and 
surface avalanches. It has been found that the models 
are related and describe the pile surface dynamics on 
different spatio-temporal scales. 

BCRE-type models may be used to simulate the fast 
processes, such as the initiation, spreading, and settling 
down of an avalanche. To describe the much slower dy- 
namics of the mean shape of a pile, the model may often 
be simplified by employing a quasistationary equation for 
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the rolling grains layer. Such a model is able to predict 
some peculiarities of small pile shapes (3^] and was re- 
cently employed for simulating the nonlinear dynamics of 
sand ripples Q . 

Unlike the BCRE models, the variational model of pile 
growth does not permit the discharged grains to roll upon 
subcritical slopes and is therefore unable to account for 
such features of sand surface as sand-ripple instability 
or surface slope deviation from the critical angle near 
the bottom of a conical pile. Indeed, these effects are 
determined by rolling of particles upon the subcritical 
slopes and are exhibited on the length scale comparable 
to the mean path of a particle prior to its being trapped. 

On the other hand, sand ripples on the dune surface 
or tiny tails at the bottom of a pile are seen only from 
a short distance. These small details are difficult to dis- 
tinguish watching from a larger distance allowing one to 
follow the evolution of a big dune or formation of a large 
pile. In such situations the BCRE model contains an- 
other small parameter. This complicates simulations and 
makes them inefficient. As has been shown in our work, 
in the long-scale limit the BCRE- type models converge to 
the variational model of pile growth. The latter model is 
more appropriate for simulating the pile surface dynam- 
ics on a large spatio-temporal scale. 



[1] L. Prigozhin, Zh. Vychisl. Mat. Mat. Fiz. 26, 1072 
(1986). 

[2] L. Prigozhin, Phys. Rev. E, 49, 1161 (1994). 
[3] L. Prigozhin, Euro. J. Appl. Math. 7, 225 (1996). 
[4] L. Prigozhin, Chem. Eng. Sci. 48, 3647 (1993). 
[5] T. Elperin and A. Vikhansky, Phys. Rev. E 55, 5785 
(1997). 

[6] H. Puhl, Physica A 182, 295 (1992). 

[7] K. Smith, Environmental Hazards (Routledge, London - 

New York, 1996). 
[8] G. Aronsson, L.C. Evans, and Y. Wu, J. Differ. Equations 

131, 304 (1996). 
[9] L.C. Evans, M. Feldman, and R.F. Gariepy, J. Differ. 

Equations 137, 166 (1997). 
[10] L.C. Evans and F. Rezakhanlou, Commun. Math. Phys. 

197, 325 (1998). 
[11] P. Bak, C. Tang, and K. Wisenfeld, Phys. Rev. Lett. 59, 

381 (1987). 

[12] J.-P. Bouchaud, M.E. Cates, J. Ravi Prakash, and S.F. 

Edwards, J. Phys. I France 4, 1383 (1994). 
[13] J.-P. Bouchaud, M.E. Cates, J. Ravi Prakash, and S.F. 

Edwards, Phys. Rev. Lett. 74, 1982 (1995). 
[14] P.-G. de Gennes, C.R. Acad. Sci., Ser. II-B, 321, 501 

(1995). 

[15] T. Boutreux and P.-G. de Gennes, C.R. Acad. Sci., Ser. 

II-B, 325, 85 (1997); 
[16] P.-G. de Gennes, Avalanches of granular materials, in: R. 

Behringer and J.T. Jenkins (eds), Powders and Grains 97 



(Balkema, Rotterdam, 1997). 
[17] T. Boutreux, E. Raphael, and P.-G. de Gennes, Phys. 

Rev. E 58, 4692 (1998); P.-G. de Gennes, Physica A 

261, 267 (1998); A. Aradian and P.-G. de Gennes, Phys. 

Rev. E 60, 2009 (1999). 
[18] L. Mahadevan and Y. Pomeau, Europhys. Letters 46, 

595 (1999); T. Emig, P. Claudine, and J.-P. Bouchaud, 

Europhys. Letters, to appear. 
[19] J.-P. Bouchaud and M.E. Cates, Granular Matter 1, 101 

(1998) . 

[20] A. Daerr and S. Douady, Nature (London) 399, 241 

(1999) . 

[21] H.M. Jaeger and SR. Nagel, Science 255, 1523 (1992). 

[22] G. Duvout and J.-L. Lions, Les Inequations en Mecanique 
et en Physique (Dunod, Paris, 1972). 

[23] R. Glowinski, J.-L. Lions, and R. Tremoliers, Anal- 
yse Numerique des Inequations Variationnelles (Dunod, 
Paris, 1976). 

[24] L. Prigozhin, Free Boundary Problems News, no. 10, 2 

(1996) . 

[25] S.N. Dorogovtsev and J.F.F. Mendes, Phys. Rev. Lett. 

83, 2946 (1999). 
[26] L. Prigozhin, Phys. Rev. E 60, 729 (1999). 
[27] K.P. Hadeler and C. Kuttler, Granular Matter 2, 9 

(1999). 

[28] F.-X. Riguidel, R. Jullien, G.H. Ristow, A. Hansen, and 
D. Bideau, J. Phys. I France 4, 261 (1994); S. Dippel, 
G.G. Batrouni, and D.E. Wolf, Phys. Rev. E 56, 3645 

(1997) ; C. Henrique et al., Phys. Rev. E 57, 4743 (1998). 
[29] M.A. Aguirre et al, Powder Technology 92, 75 (1997). 
[30] O. Pouliquen, Physics of Fluids 11, 542 (1999). 

[31] L. Prigozhin and B. Zaltzman, in preparation. 
[32] J.J. Alonso and H.J. Herrmann, Phys. Rev. Lett. 76, 
4911 (1996). 

[33] R.A. Slater, Stockpiling and Reclaiming Systems, in: 
M.C. Hawk (ed.), Bulk Materials Handling, v. 1 (Uni- 
versity of Pittsbourhg, Pittsbourhg, 1971). 



G 



